This is a toy example of application of Bayesian Optimization. We will use the Bayesian Optimization to find the maximum value of one simple function.
### install.packages("GPfit")
library(GPfit)
library(dplyr)
library(ggplot2)
library(purrr)
We will optimize this simple function:
\[\left(6x-2\right)^2sin\left(12x-4\right)\]
First we define the function in R, then draw the graph of this function in the axes, and generate n0 initial points.
### Generate the function
f <- function(x) {
return((6*x-2)^2*sin(12*x-4))
}
### Draw the graph of this function
x <- seq(0, 1, by=0.0001)
y <- f(x)
graph_data <- data.frame(x, y)
graph_data %>% ggplot() +
geom_point(mapping = aes(x = x, y = y), size = 0.01, color = "navyblue") +
theme_bw()
y_min <- min(y)
min <- 0
for (i in x) {
if(f(i) < min) {
min <- f(i)
x_min <- i
}
}
### Generate the initial points
evaluation <- data.frame("x" = c(0, 1/3, 2/3, 1), "y" = f(c(0, 1/3, 2/3, 1)))
Here we build the initial GP model. We choose the power exponential correlation function here.
fit <- GP_fit(
X = evaluation[, "x"],
Y = evaluation[, "y"],
corr = list(type = "exponential", power = 1.95)
)
x_new <- seq(0, 1, length.out = 100000)
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
mu <- pred$Y_hat
sigma <- sqrt(pred$MSE)
pred <- pred %>% data.frame()
pred %>% ggplot() +
geom_line(mapping = aes(x = complete_data.xnew.1, y = Y_hat), color = "gold") +
geom_ribbon(mapping = aes(x = complete_data.xnew.1, y = Y_hat, xmin = 0, xmax = 1, ymin = Y_hat-sqrt(MSE), ymax = Y_hat+sqrt(MSE)), fill = "navyblue", alpha = 0.3) +
geom_point(mapping = aes(x = evaluation[1, "x"], y = evaluation[1, "y"]), color = "navyblue") +
geom_point(mapping = aes(x = evaluation[2, "x"], y = evaluation[2, "y"]), color = "navyblue") +
geom_point(mapping = aes(x = evaluation[3, "x"], y = evaluation[3, "y"]), color = "navyblue") +
geom_point(mapping = aes(x = evaluation[4, "x"], y = evaluation[4, "y"]), color = "navyblue") +
xlab("x") +
ylab("f(x)") +
theme_bw()
Then we will use diferent acquisition functions to find the x to evaluate next.
### y_best is the current best y value
y_best <- min(evaluation[, "y"])
### Make the f(x) - x plot a function
plot <- function(x_next) {
pred %>% ggplot() +
geom_line(mapping = aes(x = complete_data.xnew.1, y = Y_hat), color = "gold") +
geom_ribbon(mapping = aes(x = complete_data.xnew.1, y = Y_hat, xmin = 0, xmax = 1, ymin = Y_hat-sqrt(MSE), ymax = Y_hat+sqrt(MSE)), fill = "navyblue", alpha = 0.3) +
geom_point(mapping = aes(x = evaluation[1, "x"], y = evaluation[1, "y"]), color = "navyblue") +
geom_point(mapping = aes(x = evaluation[2, "x"], y = evaluation[2, "y"]), color = "navyblue") +
geom_point(mapping = aes(x = evaluation[3, "x"], y = evaluation[3, "y"]), color = "navyblue") +
geom_point(mapping = aes(x = evaluation[4, "x"], y = evaluation[4, "y"]), color = "navyblue") +
geom_vline(mapping = aes(xintercept = x_next), color = "red", linetype = "dashed") +
geom_point(mapping = aes(x = x_next, y = Y_hat[ceiling(x_next/0.01)]), color = "red") +
xlab("x") +
ylab("f(x)") +
theme_bw()
}
1.Probability of improvement
probability_improvement <- map2_dbl(
mu,
sigma,
function(m, s) {
if (s == 0) return(0)
else {
poi <- pnorm((y_best - m) / s)
# poi <- 1 - poi (if maximizing)
return(poi)
}
}
)
x_next_POI <- x_new[which.max(probability_improvement)]
ggplot() +
geom_line(mapping = aes(x = x_new, y = probability_improvement), color = "navyblue") +
geom_vline(mapping = aes(xintercept = x_next_POI), color = "red", linetype = "dashed") +
xlab("x") +
theme_bw()
plot(x_next_POI)
expected_improvement <- map2_dbl(
mu, sigma,
function(m, s) {
if (s == 0) return(0)
gamma <- (y_best - m) / s
phi <- pnorm(gamma)
return(s * (gamma * phi + dnorm(gamma)))
}
)
x_next_EI <- x_new[which.max(expected_improvement)]
ggplot() +
geom_line(mapping = aes(x = x_new, y = expected_improvement), color = "navyblue") +
geom_vline(mapping = aes(xintercept = x_next_EI), color = "red", linetype = "dashed") +
xlab("x") +
theme_bw()
plot(x_next_EI)
kappa <- 2 # tunable
lower_confidence_bound <- mu - kappa * sigma
# if maximizing: upper_confidence_bound <- mu + kappa * sigma
### Notice here is min, not max.
x_next_LCB <- x_new[which.min(lower_confidence_bound)]
ggplot() +
geom_line(mapping = aes(x = x_new, y = lower_confidence_bound), color = "navyblue") +
geom_vline(mapping = aes(xintercept = x_next_LCB), color = "red", linetype = "dashed") +
xlab("x") +
theme_bw()
plot(x_next_LCB)
In the last section, we discussed how Bayesian optimization works in one loop. In this section, we will continue looping until we got the optimized value of the simple function.
Firstly we will choose a “good” acquisition function.
This acquisition function measures the likelihood that the next x is better(lower or higher depends on our goal) than the current optimal value.
evaluation_t <- evaluation
fit_t <- fit
t <- 1
data <- data.frame()
data$Y_hat <- c()
data$MSE <- c()
data$iteration <- c()
data$x_x <- c()
data$y_y <- c()
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
pred <- pred %>% data.frame() %>% select(c("Y_hat", "MSE", "complete_data.xnew.1"))
pred <- pred %>% mutate("iteration" = 0) %>% mutate("x_x" = NA) %>% mutate("y_y" = NA)
data <- rbind(data, pred)
acq_data = data.frame()
probability_improvement <- probability_improvement %>%
as.data.frame() %>%
mutate("iteration" = 1) %>%
mutate("x_next_POI" = NA)
colnames(probability_improvement) <- c("probability_improvement", "iteration", "x_next_POI")
acq_data <- rbind(acq_data, probability_improvement)
x_next_df_POI <- data.frame(x_next_POI) %>% mutate("probability_improvement" = NA) %>% mutate("iteration" = 1)
colnames(x_next_df_POI) <- c("x_next_POI", "probability_improvement", "iteration")
acq_data <- rbind(acq_data, x_next_df_POI)
while(t <= 10) {
evaluation[nrow(evaluation)+1,] = c(x_next_POI, f(x_next_POI))
fit <- GP_fit(
X = evaluation[, "x"],
Y = evaluation[, "y"],
corr = list(type = "exponential", power = 1.95)
)
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
mu <- pred$Y_hat
sigma <- sqrt(pred$MSE)
pred <- pred %>% data.frame() %>% select(c("Y_hat", "MSE", "complete_data.xnew.1"))
pred <- pred %>% mutate("iteration" = t) %>% mutate("x_x" = NA) %>% mutate("y_y" = NA)
data <- rbind(data, pred)
probability_improvement <- map2_dbl(
mu,
sigma,
function(m, s) {
if (s == 0) return(0)
else {
poi <- pnorm((y_best - m) / s)
# poi <- 1 - poi (if maximizing)
return(poi)
}
}
)
x_next_POI <- x_new[which.max(probability_improvement)]
probability_improvement <- probability_improvement %>% as.data.frame()
colnames(probability_improvement) <- "probability_improvement"
probability_improvement <- probability_improvement %>%
mutate("iteration" = t+1) %>%
mutate("x_next_POI" = NA)
acq_data <- rbind(acq_data, probability_improvement)
x_next_df <- data.frame(x_next_POI) %>% mutate("probability_improvement" = NA) %>% mutate("iteration" = t+1)
acq_data <- rbind(acq_data, x_next_df)
t <- t+1
}
t <- 1
nr <- nrow(data)
for (i in 4:14) {
for (j in 1:i) {
data[nr + t,] = c(NA, NA, NA, i-4, evaluation[j, 1], evaluation[j, 2])
t <- t+1
}
}
plot_1 <- data %>% ggplot() +
geom_line(mapping = aes(x = complete_data.xnew.1, y = Y_hat), color = "gold") +
geom_ribbon(mapping = aes(x = complete_data.xnew.1, y = Y_hat, xmin = 0, xmax = 1, ymin = Y_hat-sqrt(MSE), ymax = Y_hat+sqrt(MSE)), fill = "navyblue", alpha = 0.3) +
geom_point(mapping = aes(x = x_x, y = y_y)) +
facet_wrap(~iteration) +
xlab("x") +
ylab("f(x)") +
theme_bw()
x_new_df <- rep(c(x_new, NA),10) %>% as.data.frame()
colnames(x_new_df) <- "x_new"
acq_data <- cbind(acq_data %>% filter(iteration<11), x_new_df)
plot_2 <- acq_data %>% ggplot() +
geom_line(mapping = aes(x = x_new, y = probability_improvement), color = "navyblue") +
geom_vline(mapping = aes(xintercept = x_next_POI), color = "red", linetype = "dashed") +
xlab("x") +
facet_wrap(~ iteration) +
theme_bw()
print(plot_1)
print(plot_2)
This acquisition function takes how large the improvement is into account.
evaluation <- evaluation_t
fit <- fit_t
t <- 1
data <- data.frame()
data$Y_hat <- c()
data$MSE <- c()
data$iteration <- c()
data$x_x <- c()
data$y_y <- c()
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
pred <- pred %>% data.frame() %>% select(c("Y_hat", "MSE", "complete_data.xnew.1"))
pred <- pred %>% mutate("iteration" = 0) %>% mutate("x_x" = NA) %>% mutate("y_y" = NA)
data <- rbind(data, pred)
acq_data = data.frame()
expected_improvement <- expected_improvement %>%
as.data.frame() %>%
mutate("iteration" = 1) %>%
mutate("x_next_EI" = NA)
colnames(expected_improvement) <- c("expected_improvement", "iteration", "x_next_EI")
acq_data <- rbind(acq_data, expected_improvement)
x_next_df_EI <- data.frame(x_next_EI) %>% mutate("expected_improvement" = NA) %>% mutate("iteration" = 1)
colnames(x_next_df_EI) <- c("x_next_EI", "expected_improvement", "iteration")
acq_data <- rbind(acq_data, x_next_df_EI)
while(t <= 10) {
evaluation[nrow(evaluation)+1,] = c(x_next_EI, f(x_next_EI))
fit <- GP_fit(
X = evaluation[, "x"],
Y = evaluation[, "y"],
corr = list(type = "exponential", power = 1.95)
)
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
mu <- pred$Y_hat
sigma <- sqrt(pred$MSE)
pred <- pred %>% data.frame() %>% select(c("Y_hat", "MSE", "complete_data.xnew.1"))
pred <- pred %>% mutate("iteration" = t) %>% mutate("x_x" = NA) %>% mutate("y_y" = NA)
data <- rbind(data, pred)
expected_improvement <- map2_dbl(
mu, sigma,
function(m, s) {
if (s == 0) return(0)
gamma <- (y_best - m) / s
phi <- pnorm(gamma)
return(s * (gamma * phi + dnorm(gamma)))
}
)
x_next_EI <- x_new[which.max(expected_improvement)]
expected_improvement <- expected_improvement %>% as.data.frame()
colnames(expected_improvement) <- "expected_improvement"
expected_improvement <- expected_improvement %>%
mutate("iteration" = t+1) %>%
mutate("x_next_EI" = NA)
acq_data <- rbind(acq_data, expected_improvement)
x_next_df <- data.frame(x_next_EI) %>% mutate("expected_improvement" = NA) %>% mutate("iteration" = t+1)
acq_data <- rbind(acq_data, x_next_df)
t <- t+1
}
t <- 1
nr <- nrow(data)
for (i in 4:14) {
for (j in 1:i) {
data[nr + t,] = c(NA, NA, NA, i-4, evaluation[j, 1], evaluation[j, 2])
t <- t+1
}
}
plot_1 <- data %>% ggplot() +
geom_line(mapping = aes(x = complete_data.xnew.1, y = Y_hat), color = "gold") +
geom_ribbon(mapping = aes(x = complete_data.xnew.1, y = Y_hat, xmin = 0, xmax = 1, ymin = Y_hat-sqrt(MSE), ymax = Y_hat+sqrt(MSE)), fill = "navyblue", alpha = 0.3) +
geom_point(mapping = aes(x = x_x, y = y_y)) +
facet_wrap(~iteration) +
xlab("x") +
ylab("f(x)") +
theme_bw()
x_new_df <- rep(c(x_new, NA),10) %>% as.data.frame()
colnames(x_new_df) <- "x_new"
acq_data <- cbind(acq_data %>% filter(iteration < 11), x_new_df)
plot_2 <- acq_data %>% ggplot() +
geom_line(mapping = aes(x = x_new, y = expected_improvement), color = "navyblue") +
geom_vline(mapping = aes(xintercept = x_next_EI), color = "red", linetype = "dashed") +
xlab("x") +
facet_wrap(~ iteration) +
theme_bw()
print(plot_1)
print(plot_2)
### Here we save the data to use it later.
data_EI <- data
This acquisition functions balances the area where mean(x) is large and the area where sigma(x) is large. In other words, this acquisition function trades off between exploration and exploitation.
evaluation <- evaluation_t
fit <- fit_t
t <- 1
data <- data.frame()
data$Y_hat <- c()
data$MSE <- c()
data$iteration <- c()
data$x_x <- c()
data$y_y <- c()
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
pred <- pred %>% data.frame() %>% select(c("Y_hat", "MSE", "complete_data.xnew.1"))
pred <- pred %>% mutate("iteration" = 0) %>% mutate("x_x" = NA) %>% mutate("y_y" = NA)
data <- rbind(data, pred)
acq_data = data.frame()
lower_confidence_bound <- lower_confidence_bound %>%
as.data.frame() %>%
mutate("iteration" = 1) %>%
mutate("x_next_LCB" = NA)
colnames(lower_confidence_bound) <- c("lower_confidence_bound", "iteration", "x_next_LCB")
acq_data <- rbind(acq_data, lower_confidence_bound)
x_next_df_LCB <- data.frame(x_next_LCB) %>% mutate("lower_confidence_bound" = NA) %>% mutate("iteration" = 1)
colnames(x_next_df_LCB) <- c("x_next_LCB", "lower_confidence_bound", "iteration")
acq_data <- rbind(acq_data, x_next_df_LCB)
while(t <= 8) {
evaluation[nrow(evaluation)+1,] = c(x_next_LCB, f(x_next_LCB))
fit <- GP_fit(
X = evaluation[, "x"],
Y = evaluation[, "y"],
corr = list(type = "exponential", power = 1.95)
)
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
mu <- pred$Y_hat
sigma <- sqrt(pred$MSE)
pred <- pred %>% data.frame() %>% select(c("Y_hat", "MSE", "complete_data.xnew.1"))
pred <- pred %>% mutate("iteration" = t) %>% mutate("x_x" = NA) %>% mutate("y_y" = NA)
data <- rbind(data, pred)
lower_confidence_bound <- mu - kappa * sigma
x_next_LCB <- x_new[which.min(lower_confidence_bound)]
lower_confidence_bound <- lower_confidence_bound %>% as.data.frame()
colnames(lower_confidence_bound) <- "lower_confidence_bound"
lower_confidence_bound <- lower_confidence_bound %>%
mutate("iteration" = t+1) %>%
mutate("x_next_LCB" = NA)
acq_data <- rbind(acq_data, lower_confidence_bound)
x_next_df <- data.frame(x_next_LCB) %>% mutate("lower_confidence_bound" = NA) %>% mutate("iteration" = t+1)
acq_data <- rbind(acq_data, x_next_df)
t <- t+1
}
t <- 1
nr <- nrow(data)
for (i in 4:12) {
for (j in 1:i) {
data[nr + t,] = c(NA, NA, NA, i-4, evaluation[j, 1], evaluation[j, 2])
t <- t+1
}
}
plot_1 <- data %>% ggplot() +
geom_line(mapping = aes(x = complete_data.xnew.1, y = Y_hat), color = "gold") +
geom_ribbon(mapping = aes(x = complete_data.xnew.1, y = Y_hat, xmin = 0, xmax = 1, ymin = Y_hat-sqrt(MSE), ymax = Y_hat+sqrt(MSE)), fill = "navyblue", alpha = 0.3) +
geom_point(mapping = aes(x = x_x, y = y_y)) +
facet_wrap(~iteration) +
xlab("x") +
ylab("f(x)") +
theme_bw()
x_new_df <- rep(c(x_new, NA),8) %>% as.data.frame()
colnames(x_new_df) <- "x_new"
acq_data <- cbind(acq_data %>% filter(iteration < 9), x_new_df)
plot_2 <- acq_data %>% ggplot() +
geom_line(mapping = aes(x = x_new, y = lower_confidence_bound), color = "navyblue") +
geom_vline(mapping = aes(xintercept = x_next_LCB), color = "red", linetype = "dashed") +
xlab("x") +
facet_wrap(~ iteration) +
theme_bw()
print(plot_1)
print(plot_2)
Here we give a detailed example to show the Expected Improvement acqusition function.
t <- 1
while(t < 10) {
print(data_EI %>% filter(iteration == t) %>% ggplot() +
geom_line(mapping = aes(x = complete_data.xnew.1, y = Y_hat), color = "gold") +
geom_ribbon(mapping = aes(x = complete_data.xnew.1, y = Y_hat, xmin = 0, xmax = 1, ymin = Y_hat-sqrt(MSE), ymax = Y_hat+sqrt(MSE)), fill = "navyblue", alpha = 0.3) +
geom_point(mapping = aes(x = x_x, y = y_y)) +
xlab("x") +
ylab("f(x)") +
labs(paste("The", t, "iteration")) +
theme_bw())
t <- t+1
}
Then we pick the x value that has the minimal f(x) value from what we evaluated.
### Pick the x that has the minimal f(x) value from what we evaluated.
cat("The opmized x we find is", evaluation[which.min(evaluation$y),]$x)
## The opmized x we find is 0.7579476
### Compare it with the real x that minimize f(x) in [0,1]
cat("The real optimized x is", x_min)
## The real optimized x is 0.7572
We can see that we use 5 evaluations to get the x, and the result is pretty close to the real one.
In this section, we will compare the Bayesian optimization with quasi-Newton method. We will try several initial guesses.
### Firstly, we build several equally spaced initial guesses in [0, 1].
ini_guess <- seq(0, 1, length.out = 9)
### Then we will use quasi-Newton method to try to find the optimal f(x) value.
n=1
for(t in ini_guess) {
res <- optim(t, f, method = "BFGS", lower = 0, upper = 1)
writeLines(paste("The", n, "initial guess x0 = ", t, ":\nx=", res$par,"\nf(x)=", res$value, "\nTimes of iteration=", res$counts[1], "\n-----------------------------------------------"))
n <- n+1
}
## The 1 initial guess x0 = 0 :
## x= 0.142589801857248
## f(x)= -0.986325406263923
## Times of iteration= 49
## -----------------------------------------------
## The 2 initial guess x0 = 0.125 :
## x= 0.142591747065892
## f(x)= -0.986325405324459
## Times of iteration= 19
## -----------------------------------------------
## The 3 initial guess x0 = 0.25 :
## x= 0.142588440215961
## f(x)= -0.986325406235592
## Times of iteration= 10
## -----------------------------------------------
## The 4 initial guess x0 = 0.375 :
## x= 0.142591304993005
## f(x)= -0.986325405639193
## Times of iteration= 19
## -----------------------------------------------
## The 5 initial guess x0 = 0.5 :
## x= 0.142589561542382
## f(x)= -0.986325406299975
## Times of iteration= 47
## -----------------------------------------------
## The 6 initial guess x0 = 0.625 :
## x= 0.757247331284393
## f(x)= -6.02074005468035
## Times of iteration= 9
## -----------------------------------------------
## The 7 initial guess x0 = 0.75 :
## x= 0.757247332200751
## f(x)= -6.02074005468175
## Times of iteration= 6
## -----------------------------------------------
## The 8 initial guess x0 = 0.875 :
## x= 0.757247564706613
## f(x)= -6.02074005500689
## Times of iteration= 26
## -----------------------------------------------
## The 9 initial guess x0 = 1 :
## x= 0.75724733297612
## f(x)= -6.02074005468293
## Times of iteration= 10
## -----------------------------------------------
We can see there are two disadvantages of quasi-Newton method. The first one is we have to choose a “good” initial guess to find the global optimum. In our example, if the initial guess is bigger than 0.5, the quasi-Newton method will find a local optimum.
The second one is it will need more times of iteration to get the optimal value. The minimum number of iteration in our example is 13 times, which is 3 times more than our Bayesian optimization example.
Also, because here we use a simple function. If we want to optimize a complex function, which cannot or very difficult to get the derivatives, the quasi-Newton method cannot or hard to find the optimum, but the Bayesian optimization can do that.
In the previous part, we use 4 evenly spaced initial points. In this part, we will try to use different sets of initial points.
Firstly, we will try 3 evenly spaced points between 0 and 1 as the initial set of points.
evaluation <- data.frame("x" = c(0, 1/2, 1), "y" = f(c(0, 1/2, 1)))
pip <- function(evaluation, num) {
fit <- GP_fit(
X = evaluation[, "x"],
Y = evaluation[, "y"],
corr = list(type = "exponential", power = 1.95)
)
t <- 1
data <- data.frame()
data$Y_hat <- c()
data$MSE <- c()
data$iteration <- c()
data$x_x <- c()
data$y_y <- c()
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
pred <- pred %>% data.frame() %>% select(c("Y_hat", "MSE", "complete_data.xnew.1"))
pred <- pred %>% mutate("iteration" = 0) %>% mutate("x_x" = NA) %>% mutate("y_y" = NA)
data <- rbind(data, pred)
acq_data = data.frame()
expected_improvement <- expected_improvement %>%
as.data.frame() %>%
mutate("iteration" = 1) %>%
mutate("x_next_EI" = NA)
colnames(expected_improvement) <- c("expected_improvement", "iteration", "x_next_EI")
acq_data <- rbind(acq_data, expected_improvement)
x_next_df_EI <- data.frame(x_next_EI) %>% mutate("expected_improvement" = NA) %>% mutate("iteration" = 1)
colnames(x_next_df_EI) <- c("x_next_EI", "expected_improvement", "iteration")
acq_data <- rbind(acq_data, x_next_df_EI)
while(t <= 10) {
evaluation[nrow(evaluation)+1,] = c(x_next_EI, f(x_next_EI))
fit <- GP_fit(
X = evaluation[, "x"],
Y = evaluation[, "y"],
corr = list(type = "exponential", power = 1.95)
)
pred <- predict.GP(fit, xnew = data.frame(x = x_new))
mu <- pred$Y_hat
sigma <- sqrt(pred$MSE)
pred <- pred %>% data.frame() %>% select(c("Y_hat", "MSE", "complete_data.xnew.1"))
pred <- pred %>% mutate("iteration" = t) %>% mutate("x_x" = NA) %>% mutate("y_y" = NA)
data <- rbind(data, pred)
expected_improvement <- map2_dbl(
mu, sigma,
function(m, s) {
if (s == 0) return(0)
gamma <- (y_best - m) / s
phi <- pnorm(gamma)
return(s * (gamma * phi + dnorm(gamma)))
}
)
x_next_EI <- x_new[which.max(expected_improvement)]
expected_improvement <- expected_improvement %>% as.data.frame()
colnames(expected_improvement) <- "expected_improvement"
expected_improvement <- expected_improvement %>%
mutate("iteration" = t+1) %>%
mutate("x_next_EI" = NA)
acq_data <- rbind(acq_data, expected_improvement)
x_next_df <- data.frame(x_next_EI) %>% mutate("expected_improvement" = NA) %>% mutate("iteration" = t+1)
acq_data <- rbind(acq_data, x_next_df)
t <- t+1
}
t <- 1
nr <- nrow(data)
for (i in (num):(num+10)) {
for (j in 1:i) {
data[nr + t,] = c(NA, NA, NA, i-num, evaluation[j, 1], evaluation[j, 2])
t <- t+1
}
}
plot_1 <- data %>% ggplot() +
geom_line(mapping = aes(x = complete_data.xnew.1, y = Y_hat), color = "gold") +
geom_ribbon(mapping = aes(x = complete_data.xnew.1, y = Y_hat, xmin = 0, xmax = 1, ymin = Y_hat-sqrt(MSE), ymax = Y_hat+sqrt(MSE)), fill = "navyblue", alpha = 0.3) +
geom_point(mapping = aes(x = x_x, y = y_y)) +
facet_wrap(~iteration) +
xlab("x") +
ylab("f(x)") +
theme_bw()
x_new_df <- rep(c(x_new, NA),10) %>% as.data.frame()
colnames(x_new_df) <- "x_new"
acq_data <- cbind(acq_data %>% filter(iteration < 11), x_new_df)
plot_2 <- acq_data %>% ggplot() +
geom_line(mapping = aes(x = x_new, y = expected_improvement), color = "navyblue") +
geom_vline(mapping = aes(xintercept = x_next_EI), color = "red", linetype = "dashed") +
xlab("x") +
facet_wrap(~ iteration) +
theme_bw()
print(plot_1)
print(plot_2)
return(data)
}
data <- pip(evaluation, 3)
Here we use four randomly located points instead of four evenly spaced points as the initial points set.
set.seed(1234)
x <- sample.int(100,4)/100
evaluation <- data.frame("x" = x, "y" = f(x))
data <- pip(evaluation, 4)
Most of the codes are from this website!!!